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ABSTRACT 

Despite greatly improved observational methods, the presence of magnetic fields at cosmo- 
logical scales and their role in the process of large-scale structure formation still remains 
unclear. In this paper we want to address the question how the presence of a hypothetical 
primordial magnetic field on large scales influences the cosmic structure formation in numer- 
ical simulations. As a tool for carrying out such simulations, we present our new numerical 
code AMIGA. It combines an A^-body code with an Eulerian grid-based solver for the full set of 
MHD equations in order to conduct simulations of dark matter, baryons and magnetic fields in 
a self-consistent way in a fully cosmological setting. Our numerical scheme includes effective 
methodes to ensure proper capturing of shocks and highly supersonic flows and a divergence- 
free magnetic field. The high accuracy of the code is demonstrated by a number of numerical 
tests. We then present a series of cosmological MHD simulations and confirm that, in order 
to have a significant effect on the distribution of matter on large scales, the primordial mag- 
netic field strength would have to be significantly higher than the current observational and 
theoretical constraints. 
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1 INTRODUCTION 

Today, we arrived at an era where cosmology has finally reached 
the stage of a precision science. The cosmological parameters have 
been determined to a typical precision of very few percents, result- 
ing in the standard ACDM model of cosmology (Komats u et al.| 
2009), and in this context the process of the cosmic structure for- 
mation can be studied in great detail. The non-linear nature of the 
gravitational dynamics and gas physics make the problem of struc- 
ture formation virtually intractable analytically, and therefore the 
field relies on numerical simulations, which have been the driving 
force behind much of the theoretical progress. While the first codes 
were just able to follow the evolution of dark matter with the N- 
body metho d (e.g.|Klypin & Shandarin|1983||Efstathiou et al.|1985| 
|Davis et al.||1985UBarnes & Hut||1986| |Viliumsen||1989| |CouchTl 
|man|1991l|Suisalu & Saar|1995l|Kravtsov et al.|1997l|Knebe etaf] 
|2001| >, the tremendous increase in computational power in the last 
years made it possible to include more and more different com- 
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formation, AGN feedback, radiative transfer and magnetic fields in 
modern cosmological codes opens the way to study many different 
aspects of the structure formation process. 

Magnetic fields play an important role in astrophysical phe- 
nomena on many different scales. Since most of the visible matter 
in the universe is ionized and magnetic fields are found on every 
scale where they can be observed, it is natural that they could also 
play a cosmological role. Unfortunately, magnetic fields on scales 
larger than individual galaxies are much more difficult to observe 
- any measurement of magnetic fields must rely on the presence 
of radiation and a magnetized medium. So far, the largest-scale 
observable magnetic fields are inside the atmospheres of galaxy 
clusters ( |Carilli & Taylor|2002||Govoni & Feretti|2004] l, reaching 
strengths of the order of in the core regions. Detection methods 
include studies of radio synchrotron and inverse Compton X-ray 
emission from clusters ( |Harris & Grin dlay 1979, Rephael i et al.| 
|1987[ > and surveys of Faraday rotation measures of polarized ra- 
dio sources passing through the cluster atmosphere I Clar ke et al.| 
|2001[ >. However, it could be possible that there are magnetic fields 
on even larger scales than the largest observable objects. In fact, 
all the "empty space" in the universe could be magnetized l |Kron-| 
|berg et al.|f 9 99 1. A truly cosmological magnetic field would not be 
associated with collapsing or gravitationally bound structures, and 
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would be coherent on scales greater than the largest known struc- 
tures (~ 100 Mpc) or even the Hubble radius, permeating the whole 
universe. Although a new generation of highly sensitive radio tele- 
scopes like LOFAR and SKA is underway, the detection of such 
a cosmological field will probably stay out of reach for the next 
decades. Currently, it is only possible to estimate an upper limit 
for the strength of such a large-scale field (Vallee 1990; Kr onberg] 
|1994j >, which should not be higher than ~ 1 nG. In order to learn 
more about the nature and effects of magnetic fields, we have to rely 
on theoretical models and numerical simulations. Even if we may 
not be able to directly prove the existence of a large-scale magnetic 
field, the subject has important cosmological implications that must 
be considered. A large-scale magnetic field can have a significant 
impact on the dynamics of cosmic baryon flows, the thermal and 
ionization history of the universe, and the onset of structure forma- 
tion ( |Sethi et al.|2008> . 

Different theories exist on the origin of a cosmological mag- 
netic field. One class of theories suggests that the creation of a 
universal, "primordial" magnetic field happened already during a 
very early stage of the evolution of the universe. Unfortunately, at 
present such theories are highly parameter-dependent and rather in- 
conclusive ( Subramanian 2008 ), and they do not yet allow to derive 
the field strength of such a primordial field; it is currently only pos- 
sible to estimate some upper limit. Again, numerical simulations 
seem to be a very good alternative to learn more about this subject. 

Cosmological simulations including magnetic fields have al- 
ready been conducted during the last decade. It has been shown 
by numerical simulations that indeed a large-scale primordial field 
of order ~ 1 nG is needed to explain the presence of the observed 
magnetic fields in galaxy clusters (Dolag et al. 19991. There have 
also been simulations of magnetic fields in filaments (Bruggen 
|et al.|[2003|, cosmolo gical simulations studying cosmic-ray elec- 
trons jMiniati et al.||2001) , and the influence of magnetic pres- 
sure on the growth of baryonic structures l |Gazzola et al.||2007| >. 
However, aside from Dola g et al.H1999) , these earlier codes only 
included magnetic fields passively, neglecting any possible back- 
reaction effects on the baryons; or, as in Gazzola et al. (2007}, the 
magnetic component was modelled simply by adding an additional, 
isotropic pressure term in the hydrodynamic equations to account 
for the magnetic pressure. 

In order to track magnetic fields, baryons and dark matter si- 
multaneously in a self-consistent way inside a cosmological frame- 
work, it is necessary to numerically solve the full set of equations 
of cosmological magnetohydrodynamics (MHD). Codes capable of 
this task have started to be developed only very recently. They in- 
clude grid codes (Fromang et al. 2006; |LTet al.|20 08, Collin s"et al.| 
|2009] > as well as an SPH code (jpolag & Stasysz yn|2009[ l, none of 
which were publicly available at the time of writing. In this paper, 
we present the new cosmological MHD code AMIGA, aimed to close 
this gapQ 

The AMIGA code originally started as a pure ,/V-body code to 
study dark matter structure formation. It is the successor of the 
MLAPM code, a very powerful and memory-efficient AMR code pub- 
lished by Knebe et al. (200TJ. Here, we present a new numerical 
solver for cosmological MHD, now implemented into the AMIGA 
code. It greatly improves the possibilities of the code, allowing to 
model dark matter, baryons and magnetic fields simultaneously in a 
fully cosmological setting. The code utilizes the transformation to 
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supercomoving coordinates, which greatly simplifies the numerical 
solution of cosmological MHD equations. There are implemented 
techniques to properly resolve strong Shockwaves and supersonic 
flows in the baryon component, and to ensure the important condi- 
tion of a divergence-free magnetic field down to machine precision. 
We also present a series of test problems, in order to verify the high 
accuracy of the code. 

After a technical description of the underlying principles and 
numerical methods, we use this new powerful tool to investigate 
and quantify the influence of a primordial magnetic field on the 
cosmic structure formation on large scales. Recent numerical ef- 
forts in cosmological MHD have concentrated on the modelling of 
magnetic fields inside individual galaxy clusters, since they are di- 
rectly observable. For example, Dubois & Teyssier (2008) focused 
specifically on the magnetic field inside one simulated cluster and 
found a relation between the field strength in the cluster core and 
cooling processes of the intracluster gas. We want to go a different 
path and study the influence of a hypothetical universal magnetic 
field, filling the whole universe, on large-scale structure formation. 
We currently only have some constraints on the maximum value 
of such a field, and no working general theory describing its ori- 
gin or its development until the present time. At this juncture, it 
seems reasonable to choose a more pragmatic strategy. We want to 
address an important question for future cosmological simulations: 
suppose a cosmological primordial field exists - could it have a dy- 
namically significant influence on the other constituents, dark mat- 
ter and baryons? Does it need to be considered when performing 
simulations of the large-scale structure? At what strengths of a pri- 
mordial field do its effects become relevant, and how do these field 
strengths compare to the current constraints? A series of numerical 
simulations of the large-scale structure formation, conducted with 
the new cosmological MHD code AMIGA and including primordial 
magnetic fields of different strength, is presented here to address 
these questions. 

The outline of this paper is as follows. Section|2]is dedicated 
to our new cosmological MHD code AMIGA. We present the super- 
comoving framework, in which we formulate the equations of ideal 
MHD {2.1) , the numerical scheme implemented in AMIGA {2.2) , and 
then we carry out different numerical tests to ensure that the code 
is functioning accurately < |2.3| >- Section [3] presents our simulations 
of structure formation with a primordial large-scale magnetic field. 
We first use the MHD pancake formation as a toy model to estimate 
the effect of the fields {}.!) , and then present our 3D cosmological 
simulations with magnetic fields and analyse the obtained numer- 
ical data {$.2) . Our conclusions are summarized in section]?] The 
derivation of the supercomoving MHD equations is given in the 
Appendix. 

2 AMIGA 

AMIGA is a cosmological grid code containing the ,/V-body solver 
with adaptive mesh refinement from its predecessor, the MLAPM 
code (Knebe et al. 2001 1, which is used for the dark matter and 
gravity equations, and a newly developed MHD solver to track the 
baryon physics and magnetic fields on a regular grid. 

2.1 Supercomoving ideal MHD equations 

Our simulations contain dark matter particles, treated by an A'-body 
code, and a baryon component that behaves like an ideal, supercon- 
ducting plasma, together with a magnetic field. To treat all these 
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components simultaneously in a self-consistent way, the equations 
of ideal magnetohydrodynamics (MHD) have to be solved together 
with the dark matter particle equations in a fully cosmological set- 
ting. For this, the equations of MHD are usually transformed to the 
comoving frame, defined by 

r 

x = - 
a 

In this frame, in the absence of additional forces, the mass 
points are at rest and the local density remains constant. However, 
if applied to the equations of MHD, the transformation renders 
them into equations with lots of additional factors explicitly de- 
pending on the cosmological expansion factor a(t); these equations 
no longer have the form of hyperbolic conservation laws. Neverthe- 
less, most other cosmological MHD codes use this formulation (Li 
|et al.|2008] |Collins et al.|2009| >. A different transformation to so- 
called supercomoving coordinates has been proposed by Martel & 
Shapiro jl998f to cast the equations into a more convenient form|j 



It is defined in a very similar way, but additionally, the physical 
time t is replaced by a new function t x depending on the expansion: 

r dt 

x = -; dt x =- (1) 
a a z 

All time derivatives are now formulated in respect to that new func- 
tion, and the equations are transformed accordingly. Here, we apply 
this transformation to the full set of MHD equations. Additionally, 
the physical quantities therein get substituted by a set of new 'su- 
percomoving' quantities: 



Px 



- pa 

- a 2 {(p + -adx 2 ) 
-a s p 



T x = a z T 
S x = a*>-*>S 



(2) 



where p is the baryon density, <p the total gravitational poten- 
tial, p the thermal baryonic pressure, e the thermal baryonic energy, 
T the temperature, S the modified entropy (definition see section 
|2.2.4| (, *H the supercomoving Hubble constant and B the magnetic 
field strength. With the supercomoving framework defined in that 
way, the substitution causes most of the a(t) depending terms to 
cancel out and results in the following equations (the x subscripts 
are dropped from here on): 
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The derivation is presented in the appendix. Equations ([3} and 
(|4) are the equations of motion for the collisionless dark matter 
(DM) particles, where x DM is the position and v DM the velocity, 
respectively. {5} is Poisson's equation for the total gravitational po- 
tential (f>, where p m is the total density of combined gas and dark 
matter, and p m is the average total density of the simulated box, 
given by 



Ptot - &0 Peril - Qo 
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Next there are the supercomoving ideal MHD equations in conser- 
vative form, where equation l|6j is the conservation law for gas den- 
sity p, equation |7| for the gas flow momentum pv, and |8| for the 
total energy density pE of the gas. B is the supercomoving mag- 
netic field, whose evolution is given by the law of induction j5|, 
subject to the divergence-free condition ( [10) . The thermal pressure 
p is obtained via an ideal equation of state, 



p = (r - i)p« 



(12) 



where the adiabatic index equals y = 5/3 for a non-relativistic, 
monoatomic ideal gas, while the internal energy density ps of the 
gas follows from the total, kinetic and magnetic energy densities: 



1 , B 2 



(13) 



Note that formally the equations |3]( to |9]( closely resemble 
their non-comoving counterparts. The only differences are in Pois- 
son's equation for the gravitational potential and the two additonal 
magnetic Hubble terms at the right side of equations {8} and j9|. 
These are now the only places where cosmology explicitly enters, 
namely in the form of the supercomoving a(t) function, which has 
to be determined depending on the adopted cosmological model. 
These properties of the supercomoving MHD equations make them 
easier to implement than comoving MHD, while still containing the 
same physics. In particular they make it very easy to employ numer- 
ical schemes originally designed for non-cosmological purposes. 

The code uses the following internal units: The distance unit is 
the comoving boxsize B , so that i,j/,ze [0, 1] always; the density 
unit is the average density < | 1 1\ , so internally 6 = p — 1; the unit 
for supercomoving time is the Hubble time I /Ho; and the magnetic 
field unit is defined by setting the magnetic constant to unity: p = 1, 
so it disappears from all equations. The expansion factor a(t) is 
evaluated by numerically integrating 



da 
df 



n A (a 2 - 1) + — - £l m + 1 
a 



1/2 



(14) 



internally in the code (this relation results from the Friedmann 
equation). Note that here, t is the supercomoving time (hence the 
additional a 2 ) and H = 1 due to the internal units. 



We like to note that codes such as, for instance, RAMSES |Teyssier|2002) 
and ART ( Kravtsov et al. 1997 2002 1 also implement supercomoving coor- 
dinates. However, neither of the code description papers has yet shown the 
full set of the corrresponding equations and their derivation as presented 
here and in the Appendix, respectively. 



2.2 Numerical scheme 

For an elaborate description of the AMR solver for equations 
(|5j we refer the reader to the MLAPM paper by |Knebe~et 



Below, we present the new solver for the MHD equations 
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2.2.1 MHD solver 

The MHD solver of AMIGA serves to solve the cosmological MHD 
equations |6]( - It consists of a second-order unsplit Godunov- 
type central scheme and a constrained transport scheme to ensure 
a divergence-free magnetic field down to machine precision. It is 
essentially an expanded, cosmological version of the solver used 
by the Nirvana code (Ziegler 2004, Ziegler 2005), which in turn 
adopts the KNP solver for hyperbolic conservation laws ( [Ku rganov 
& Petrova|20 01 )01n the following section, we present the numeri- 
cal algorithm, including the KNP solver, as implemented in AMIGA; 
for more on the theory behind the scheme, we refer the reader to 
these articles. 

The three hydrodynamic conservation laws |6](, j7| and ([8j can 
be written in general vector form: 



du + df + df + df _ 

dt dx dy dz 

where u is a vector containing the hydrodynamic variables, 
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f,f,f are the flux functions 
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and S u are the source terms 
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AMIGA was developed from a particle-mesh code and inher- 
ited its grid structure. We use the cells of this grid to locally store 
discrete values of the MHD quantities u, B. 

The hydrodynamic al quantities u are stored as cell-averaged 
values Uj j k at the centres of the grid cells i,j,k whereas the mag- 
netic field B is instead arranged in a "staggered grid", i.e. it is stored 
on the cell faces with a staggered collocation of the components 
B x , B y , B z (see Figure[T|a). Thus, every component of the B field is 
stored at another interface of the cell. 



3 We like to note that the KNP flux is equivalent to the HLL flux formula in- 
troduced by Harten et al. ( 1983 1 as two-speed approximate Riemann solver. 
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Figure 1. a) The variables stored at every cell i, j, k of the 3D grid. The 
cell-averaged hydrodynamic quantities Uuji are stored at cell centres (', j, k, 
while the components of B are stored on different cell faces ("staggered 
grid"), b) The reconstructed MHD quantities on all six cell faces. 



In three dimensions, the KNP solver requires the reconstruc- 
tion of all these quantities to all six faces of every cell (see Figure^ 
b). We will denote the six interfaces with the letters W, E,S,N, T, B. 
For example, we obtain a hydrodynamical variable u at the inter- 
faces lying in x direction (E and W interfaces) by 



U i,j.k ~ "i,j,k + ~{8xU)ijk 



ll i.j,k = "ij-k ~ ~($.x u )ijk 



(18) 
(19) 



where (S x u) is a TVD slope limiter. For cosmological MHD simu- 
lations, we use the slope limiter of | van Leer| ( |1977^ : 



(S x u)i 



For the y direction (N and S interfaces) and the z direction (T and 
B interfaces) the formulae are analogue. 

Since there are magnetic terms present in the hydrodynamic 
flux functions, we also reconstruct all magnetic field components 
at these interfaces. The only difference is that, due to the staggered 
grid collocation of B, we have to average over pairs of opposing 
cell interfaces on the way. For the interfaces lying in x direction 
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this means 
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(20) 



B z.,M =o( S c,.;,* + | +» ! i lil 4-(«)^.-M) W .) 



D ^i,j,t ~ D xi-{,j,k 



p- 1 " _ p x F' J - F' J 
d i+ij,* i-ij,* ij+j^ >'-j-* ,„„ 
My* = ; : — (26) 



At 



Ax 
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Note that no time discretization has been specified yet. We will 
later use (26) to update Uj j k applying a second-order Runge-Kutta 
scheme for the time integration (see section |2.2.3} . 



Note that the B x component does not get reconstructed, since 
it is already stored at the needed interface. Again, for the other two 
directions there are analogue expressions. 

Now, we calculate the flux functions / at E,W,N,S,T,B loca- 
tions by putting the interface values of u and B into the definition 
(16}. Once we have them, the numerical fluxes in and out of each 
cell at each interface are calculated utilizing the KNP flux formula: 



F x 
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i+ k ,j,k i+ i J,k 
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ij+ j ,k 1 ,k 
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i,iM{ Ki,k+{ 

In these flux formulae, a* denotes the maximum (+) and min- 
imum (-) local speed of the hydro density flow at the cell surface 
in x-direction (wavespeed estimate of Davi s|1988| >: 

<Uk = max{fe + c f$ f +i i j,k> (Pi + Cftj* 0) (22) 



and fe*, c* the same for y and z, respectively. The expression 



(23) 



is an upper limit for the possible characteristic wave speed in the 
medium (fast magnetosonic speed), where 



is the sound speed in the medium and 



(24) 



(25) 



is the Alfven speed. All these quantities get calculated on-the- 
fly using the reconstructed MHD variables and the equations §\2\ 
and (13). 

By adding the fluxes through all six cell interfaces, we now 
have the total flux in and out of the cell: 



2.2.2 Constrained transport ( CT) 

In order to track the time evolution of the magnetic field B as well, 
we want to solve the induction equation (9) in a similar way. But 
when introducing magnetic fields to such grid algorithms, one is 
immediately faced with the problem that the solution for B must 
comply at all times to the additional condition V ■ B = down 
to the highest possible precision. Otherwise magnetic "sources" 
(monopoles) would be introduced that would lead to unphysical 
results (like forces parallel to the field direction). Physically, V • B 
is a conserved quantity. But this is not the case for numerical cal- 
culations - a nonzero divergence will inevitably build up due to 
numerical errors, even if B was divergence-free at the beginning of 
the simulation. Even worse, numerical nonzero V • B usually grows 
exponentially (Brackbill & Barnes 1980), and the code will crash. 

There are a handful of techniques to remedy the situation 
(see Toth 2000 for a review and comparison study). Brackbill & 
|Barnes " ]l980) introduced the "divergence cleaning" (or "Hodge 
Projection") approach, which solves an extra Poisson's equation 
to recover V ■ B = at each time step. But it was found later 
that the divergence cleaning can introduce substantial amounts of 
additional spurious structur e ([Balsara & Kim|2004[ >. The second 
method ( |Powell et al.|1999|[Dedner et al.|2002) extends the MHD 
equations to produce an additional divergence wave, which then 
advects the divergence out of the domain. This generally works; 
however, in cosmological simulations we always work with peri- 
odic boundary conditions. Thus, a wave cannot leave the domain, 
and this method is not applicable. 

In AMIGA, we use the arguably most elegant solution, the con- 
strained transport (CT) method by Evans & Hawley ( 1988 1. In this 
method, the components of B are arranged in a way that ensures 
V • B = by definition. This is the reason why we introduced the 
staggered grid. 

Another issue with incorporating the induction law in the cho- 
sen finite-volume scheme is that the conserved quantity, i.e. the 
magnetic flux, is defined on a surface rather than on a volume (such 
as density). This leads to the fact that, as opposed to equations (6)- 
(8), the induction equation (9), albeit it is a conservation law, con- 
tains the curl operator Vx instead of the divergence operator V- . 
Here, it is possible to apply a trick that still allows to use the same 
numerical scheme for the magnetic field as for the hydro variables 
u. If we write E = —v x B, then the induction equation becomes 



8B 



+ V x £ = S B 



(27) 



with the magnetic Hubble source term Sb = j'HB (in the non- 
cosmological case it is zero). Now this equation can be transformed 
into divergence form using an antisymmetric flux tensor, i.e. 
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8B 



E- -E y ' 
+ V -E- -E x = S B (28) 

v Ey ~E X 

which is formally analogous to l |15[ >; but instead of flux functions 
l [16) , we use the components of the antisymmetric flux tensor. It 
is essentially just a "resorting" of the vector components to make 
the curl appear formally as a divergence. With this, it is possible to 
construct numerical fluxes E using the same KNP flux formula as 
before: 
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Now, to ensure V • B = 0, the constrained transport method 
enters. The idea is to discretize the magnetic field and the magnetic 
fluxes in such a way that the V • B = condition follows by def- 
inition of the scheme and is therefore conserved down to machine 
precision. 

From the face-centered fluxes {29} we calculate edge-centered 
fluxes. An easy way to calculate £ on a cell edge is averaging over 
the four interfaces touching this edge (see Figure[2}: 
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(30) 



The other nine edges bordering cell i, j, k are obtained in the 
same way with the according indices. It has been pointed out re- 
cently by |Gardiner & Stone| ( |2005| |2008} that this is actually not 
the best way to construct edge-centered fluxes and that for certain 
cases, a reconstruction algorithm with proper upwinding gives bet- 
ter results than the simple averaging. However, for the tests and 
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Figure 2. Face-centered and edge-centered E fields in the constrained trans- 
port scheme. 



simulations presented here, this has no relevant effects. If it be- 
comes necessary in the future to improve the algorithm, another 
scheme like the one by Gardiner & Stone (2005, 2008) may be 
implemented by simply adjusting equation (|30|> accordingly. Al- 
ternatively, a numerical dissipation term can be introduced in the 
induction equation to smear out any possible numerical noise. 

Using the edge-centered fluxes, we get the temporal change of 
the staggered magnetic field components, in analogy to the hydro 
flux (f26l: 
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(31) 



Ay 

By writing out V ■ B with these discretizations one can im- 
mediately see that d(V ■ B)/dt = by definition. Therefore, with 
compatible initial conditions, V • B = is conserved at all times. 



2.2.3 Time integration 

In order to calculate the temporal changes of the MHD variables, 
we discretize the equations {26} and {3T} in time with a standard 
second-order accurate two-step Runge-Kutta method. At a given 
time r, we have u' and B' stored as described before. First, we 
estimate the appropriate timestep At by using the usual timestep 
criteria. The time step At should not exceed the actual age of the 
universe, 



At ^ 



1 



(32) 
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the fastest dark matter particle in the box should not travel farther 
than some fraction e, of one grid cell during one timestep, 

€\ Ax 



At < 



VDM.il 



(33) 



and the same must hold for the fastest baryon flow speed encoun- 
tered in the medium, 



At < 



e 2 Ax 



(34) 



where e 2 , the so-called CFL number, should always be < 0.5 (CFL 
criterion). Now, denoting the right-hand side of equation ( |26} as 
F(u,B) and the right-hand side of equation |3TJ as E(u,B), we 
perform a predictor timestep 



«' +A " 

B , + A„ 



u' + At-F(u',B t ) 
B' + At •£(«', B') 



(35) 



These predictor values are used to calculate timestep-centered val- 
ues 

1 



B-T- = _(„' +«' +A ") 

2 



(36) 



B' 



(B' + B' +A ") 



and in a final step, u and B get stepped forward in time using these 
timestep-centered values: 

*• At , a ;+^« 

+ y • F(u ,B ) + At ■ S * 2 (37) 



B 



• E(u'- 



,B' +& ") + At-S B 



where S„ are the timestep-averaged hydro source term s l|17[ l and 
Sb = HB/2 is the magnetic Hubble term from equation {9JL 

It is important to point out that we must use a different time 
integration scheme here than the one in the Af-body part. While 
leapfrog-based integrators like the one used by the Af-body solver 
are well suited for Hamiltonian-type equations of motion and in 
particular the W-body problem, they are unstable for hyperbolic 
conservation laws like the MHD equations. However, the time in- 
tegration schemes in the Af-body solver and the MHD solver are 
connected only through the gravitational potential 0. For both time 
integrators, the timestep-averaged gravitational potential f+ T is 
needed, so we can compute that from the time-averaged total den- 

<+- n 

sity p m 2 (the sum of baryon and dark matter density). Figure 13] 
shows a flowchart of a full AMIGA timestep, illustrating how the 
A'-body-solver and the MHD algorithm are interconnected. 

The code is completely modular, i.e. it is possible to run a pure 
A'-body simulation, a pure hydrodynamic simulation, an MHD sim- 
ulation or a combination of everything. For non-cosmological runs 
like the test cases presented before, it is possible to set a = 1, a = 0, 
and the supercomoving MHD equations reduce to the ordinary 
MHD equations in proper physical coordinates. The gravity solver 
and the periodic boundary conditions can also be modified or dis- 
abled. 



2.2.4 Supersonic flows and the dual energy formalism 

When including gas physics in cosmological simulations, due to 
the extreme gravitational forces the gas flows can be accelerated 
to highly supersonic speeds, reaching Mach numbers of 100 and 
more. While the shocks and discontinuities that are created by such 
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Figure 3. Flowchart of AMIGA'S timestepping on a regular grid. The stored 
values are dark matter positions x and velocities v, the hydrodynamical 
quantities u (cell-centered values) and the magnetic field B (staggered face- 
centered values). 
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flows can be captured very well by the KNP solver, they are also 
followed by highly supersonic bulk flows of cold gas. A serious 
numerical problem occurs when trying to describe such flows with 
ideal MHD equations. 

At different places in the code, the thermal energy density 
pe and pressure p have to be calculated. Normally, this happens 
through equations \12\ and < | 1 3[ >. The problem lies in the fact that 
in such cold, highly supersonic bulk flows, the value of pe will be- 
come several orders of magnitude smaller than the total energy den- 
sity pE. Expression \13\ then contains a small difference of large 
numbers, leading to wrong results due to limited floating point pre- 
cision. The thermal pressure and therefore the gas temperature can- 
not be tracked accurately anymore. 

To remedy this situation, |Ryu et aL| ([T993 1 proposed to intro- 
duce the modified entropy as an additional equation from which 
the thermal energy could be calculated. Alternatively, Bryan et al. 
( 1995) use an equation for the thermal energy itself (which is a bit 



problematic because it is not a conservation law). We follow the 
|Ryu et aT] (l993 ) method and define the modified entropy as 



S = 



P 



(38) 



In supercomoving coordinates, the supercomoving modified en- 
tropy follows the conservation law 



dS_ 
~8t 



+ V -(Sv) = -f{S(3y-5) 



(39) 



(see appendix for the derivation), where the right-hand side equals 
zero for y = 5/3. 

For the whole simulation, we solve this equation alongside 
equations l|6j through ([8j with the KNP solver. Now, whenever the 
thermal energy cannot be calculated accurately from the usual set 



of equations 
of equations 



and i gl- 
and 1391 



the "E system" 
- the "S system" 



, we use the set 
instead. In this 



system, the pressure and thermal energy are calculated as follows 

Spt- 1 



Spf 



7-1 



(40) 



After each timestep, the two systems have to be resynchronized: 
if the S system was used, the total energy has to be updated to be 
consistent with the new internal energy; if the E system was used, 
the entropy has to be updated according to equation l |38[ >. 

The crucial step here is how to determine when to use the en- 
tropy for the calculation. A possible choice is to do this whenever 
the ratio of pe to pE gets smaller than some threshold parameter 
(e.g. in |Collins et al.|2009) : 

pE-pv 2 /2-B 2 /2 



pE 



(41) 



This works very well for most cases. However, in cosmological 
simulations sometimes another situation occurs when all energy 
components are near zero numerically, for example in the low- 
density regions between the shocks in the double pancake test (see 
section|233J. The condition is false, nevertheless the pressure 
can not be tracked accurately. In order to deal with this issue, we 
propose a new approach: Instead of using the S system only in cases 
where l |41[ > is true, we reverse the original condition and use the S 
system always, given that the entropy is conserved (that is, outside 
of shocks). Whether we are in a shock or not gets estimated by an 
additional criterion, checking for steep pressure gradients: 



|Vp| 



< m 



(42) 



Whenever either l |41[ > or \A2) is true in a grid cell, we calculate the 
thermal energy using the S system. This new method gives accu- 
rate thermal quantities not only in strong shocks, but also in very 
cold low-density regions. Of course, since technically the S system 
abandons strict energy conservation in favor of accurately tracking 
the temperature, we must make sure that the use of the S system 
does not have a dynamical effect on the other hydrodynamic vari- 
ables by choosing the parameters low enough. For the cosmological 
MHD simulations presented here, we used r)\ = 0.001 and 772 = 0.3. 

The cell-averaged value of the magnetic energy density B 2 /2 
is required here for compatibility with the other energy terms. It is 
calculated by averaging over the opposing pairs of face-centered B 
components that enclose the cell: 



1 1 



i,j,k 



(43) 



+ B 



yi,j-k,k> + ( B zi,jJc+i 



+ B 



2.3 Code testing 



To verify that the code is working correctly we applied it to a set of 
standard test problems. The iV-body solver of AMIGA comes from 
its predecessor MLAPM ( Knebe et al. 2001 1. It has been thoroughly 
tested therein and successfully used for cosmological dark matter 
simulations (e.g.|Gill et al.|2004|?||2005l|Warnick & Knebe|2006j 
|Warnick et al.|2008| >. Therefore we can concentrate here on test- 
ing the newly implemented MHD solver and its interplay with the 
gravity solver. 

The hydrodynamic part of the solver is applied on a ID test, 
the Sod shock tub e <|Sod||1978| >, and a 3D test, the Sedov-Taylor 
blast wave jSedov||1959^ . To verify the MHD solver and the CT 
scheme we use the Brio-Wu problem ([Brio & Wu||1988} and the 
Orszag-Tang vortex (Orszag & Tang 1979|. Then, combining MHD 
with the gravity solver and the cosmological expansion, we present 
the double pancake test of |Bryan et al.|(T995) , which also serves as 
a stringent test on the dual energy algorithm. 

The computational domain for all tests is x,y,z £ [0,1], con- 
forming with the internal code units. All numerical runs up to 
N = 256 cells of box length have been performed on the full three- 
dimensional N 3 box, even for ID test problems, to test the code 
under more realistic circumstances. For higher resolutions we used 
a reduced ID box to save computing time. It turned out that both 
recover the exact same result. Furthermore, in the case of pure hy- 
drodynamic tests with no magnetic field, full MHD runs with the 
initial B set to zero recover the exact same result as purely hydro- 
dynamic runs. 



2.3.1 Sod Shock tube 

For the Sod shock tube test, the simulation box is divided in two 
halves with constant initial states separated by a barrier between 
them that is removed at t = 0. This generates a strong shock wave 
moving to the lower density region, a sound wave (rarefaction) in 
the opposite direction and a contact discontinuity. This simultane- 
ous presence of different phases makes the shock tube an excellent 
method to check how well a code handles strong shocks. We chose 
the same initial conditions as in the original paper of Sod ( [1978| (. 
The left and right initial states at t = are: 



Pl ■ 
Pl 



1 



p R = 0.125 
Pr = 0.1 
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The initial velocity is zero everywhere, the polytropic index is y = 
1.4. The boundary conditions are non-periodic and the system is 
evolved until t mr i = 0.2. By then the main shock will be located 
at x = 0.85. Figure [4] shows the results perpendicular to the shock 
plane. 

In the code comparison suite of Tasker et al. ( 20081, this test is 
applied to other astrophysical codes. In direct comparison, AMIGA 
handles the problem very well. The main shock is between three 
and four cells wide, an accuracy comparable to PPM grid codes. 
All features of the analytical solution are recovered very accurately 
without oscillations or other numerical artefacts, except for a slight 
overshoot in the internal energy at the contact discontinuity. This is 
a common feature in grid codes and quickly disappears with higher 
resolution. 

The analytical reference solution for this problem was gener- 
ated with an exact Riemann solver algorithm based on Toro ( 1999 1, 
using a resolution of N = 10 4 . 



2.3.2 Sedov Blast wave 

The Sedov blast wave test is performed by injecting a large amount 
of thermal energy in a small, point-like region with uniform cold 
gas around it. This causes a strong explosion with a spherical shock 
front propagating outwards. The test is particularly useful to check 
if spherical symmetry is preserved by the code. It is important that 
on a cubic grid, shock fronts that are aligned parallel to the grid 
are resolved the same way as those moving at an oblique angle, be- 
cause otherwise we would introduce an artificial anisotropy. Also, 
the shock front is very narrow and thus numerically challenging 
to resolve. As a reference we use the known self-similar analytical 
solution jSedov|1959t . 

For this test, the gas is at rest with p = 1 and v = everywhere. 
We inject the energy E = 1 in a spherical region of radius 3.5 cells 
in the centre of the simulation box. Then, the initial pressure equals 



3(y- Wo 
4nr 3 

io- 5 



if r < 3.5 Ax 



else 



We use y = 7/5 and evolve the blast wave until / = 0.0508. 
The shock is then located at r = 0.314. According to Tasker et al. 
(2008 ), the r = 3.5 cells sphere is a good approximation of a point- 
like energy source, if we use a uniform grid with N = 256 or higher; 
so we use exactly this resolution. The results are compared with the 
analytical solution in figure|5] 

The code conserves spherical symmetry and recovers the an- 
alytical solution well. The shock front is smoothed over a width 
of approx. four grid cells, so there is not much broadening due to 
a shock propagation on different angles with respect to the grid. 
The anisotropic scatter is not larger than one grid cell. The peak 
amplitude is somewhat lower than the analytical solution, but still 
very well compared with other codes (Ricker et al. 2000; Tasker 
|et al.| [2008l. The lowering is partly due to the fact that we use a 
finite spherical region instead of a really point-like source, which 
is just impossible with a grid code. Some codes also suffer from 
other problems: the shock position is sometimes underestimated by 
as much as 4%, for example Enzo (Zeus) in the Tasker et al. ( 2008 1 
code comparison, since the initial energy lies in a region made of 
cubical cells and is therefore not exactly spherically symmetric. It 
can result in a deformed Shockwave lagging behind the analytical 
solution, and the position will be wrong. However, AMIGA does not 
suffer from such problems, even if no technique is applied to make 



the start region more spherical (e.g. Gaussian smoothing or some 
other weighted distribution), and always recovers the correct shock 
front position. 

For both hydrodynamic tests in general, we find that the nu- 
merical accuracy of AMIGA'S hydrodynamic shock capturing is on 
par with the most popular astrophysical grid codes used today. 



2.3.3 Brio-Wu problem 

Now we want to test whether the MHD equations are implemented 
correctly into the solver. We use the test of Brio & Wu (19881, 
which is one-dimensional, so the constrained transport reduces to a 
simple advection (we will move on to a multi-dimensional test af- 
terwards). The Brio-Wu test is very similar to the Sod shock tube, 
with left and right initial states and a Riemann discontinuity in be- 
tween; but in addition it features a magnetic field that has com- 
ponents both parallel and perpendicular to the shock plane, inter- 
acting with the shock and the different discontinuities. We use it 
to check the correct implementation of MHD equations and MHD 
shock capturing in one dimension, before continuing with multi- 
dimensional and cosmological tests. The initial conditions for this 
test are: 



Pi 

PL 



0.75 
1 




PR =0.125 
Pa = 0.1 
(Q.15\ 



B R 



-1 




This leads to the propagation of all seven MHD waves (2 shocks, 2 
Alfven waves, 2 slow magnetosonic waves and a contact disconti- 
nuity) to travel through the box. For these initial conditions, two of 
the waves will have almost the same speed and interfere with each 
other, causing the overshoots typical for this test. 

Again, the velocity is zero everywhere, but this time we use 
7 = 2. The system is evolved until t = 0.1. Since there is no analyt- 
ical solution known for this problem, we use a high-resolution run 
with N = 1024 as a reference. For comparison, this test can also be 
found e.g. in |Ryu &~Jo nes ( 19951. The numerical results obtained 
with AMIGA are shown in figure[6] In the pre-shock region, there is 
an additional overshoot in the x-direction velocity that disappears 
only at higher resolutions; otherwise, the results are quite accurate 
and show a good convergence with higher resolution. We acknowl- 
edge that the MHD equations are implemented correctly and the 
shock capturing is accurately handled by the code in the full MHD 
case. 



2.3.4 Orszag-Tang vortex 

The next task is to check the multidimensional MHD behaviour: the 
correct implementation of the constrained transport algorithm and 
the conservation of V • B = 0. One of the most popular benchmark 
tests for that purpose is the Orszag-Tang vortex (Orszag & Tang| 
|1979| >. This 2D test features an initially smooth flow that quickly 
develops MHD shocks and shock-shock interactions, and eventu- 
ally breaks down into supersonic MHD turbulence. The initial con- 
ditions for this test are: 
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Figure 5. Numerical solution of the Sedov blast wave. The dots correspond to cell values of the AMIGA run at t = 0.508 with a resolution of 256 3 ; the solid 
line is the reference solution computed from the analytical formulae given in Sedov 1 1959 1. 



Cosmological Magnetohydrodynamics with AMIGA 1 1 



i 

0.9 
0.8 
0.7 
p0.6 
I 0.5 
0.4 
0.3 
0.2 
0.1 



N=64 
N=256 
N=1024 



0.2 0.4 0.6 0.8 




0.2 0.4 0.6 0.8 



0.2 0.4 0.6 0.8 




Figure 6. Numerical solution of the Brio-Wu problem. Since no analytical solution is known, a high-resolution run with N = 1024 serves as the reference 
solution. This test can be found in e.g. Ryu & Jones 1 1995 1. There is an overshoot present in the x-direction velocity that disappears only at higher resolutions; 
otherwise, the results compare and converge extremely well. 
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Figure 7. Numerical solution of the Orszag-Tang vortex at t = 0.5 with N = 256 resolution. From left to right: temperature and magnetic energy density 
distribution in the x-y plane; gas pressure along a cut at y = 0.4277. 
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(44) 



We use periodic boundary conditions and y = 5/3, leading to 
c s = ^jjplp = 1 everywhere. Note that here, as for any MHD test, 



the initial magnetic field must be chosen so that it is divergence- 
free. The system is then evolved until t = 0.5 . 

The numerical results are shown in figure|7] maps of the tem- 
perature and magnetic field energy density in the computational 
plane, and the gas pressure along a cut at y = 0.4277. AMIGA recov- 
ers the characteristic, complex shape of the solution in great detail, 
including the thin thread-like structure in the middle of the box. 
Magnetic field components are tracked correctly, and most impor- 
tantly, V • B equals machine zero at all times and positions. For 
comparison, the same test performed on other MHD codes can be 
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Figure 8. The double pancake problem at z = with different resolutions (in supercomoving coordinates and code units; velocity in km/s). For this test we set 
B = and use the dual energy formalism as described in subsection |2.2.4| The boxsize is 64 Mpc/h. 
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Figure 9. Convergence check for the double pancake problem. The left plot is the relative error in the linear regime (at redshift z = 20). The right plot is the 
self-convergence error (definition see text) for the final solution at z = 0. Since the evolution is highly non-linear at this time and features singularities, the 
convergence rate is quite low at poorer resolutions (around N~°' S ). 



found e.g. in Ryu et aL|(1998[ l, |Fromang et al. (20061 (grid codes) bation with a certain wavelength A = 2n/k, it can be seen as a 
and |B0rve et al. 1 2006 i, |Dol ag & Stasyszyn ( 2009 ) (SPH codes). single-mode analysis of actual cosmological structure formation. 



2.3.5 Double pancake test 



The cosmological pancake formation (Zel'dovic h|1970| > is a very 
popular test for cosmological hydrocodes, because it combines all 
of the essential physics (hydrodynamics, cosmological expansion 
and gravity) and is a very stringent test due to the strong shocks 
and non-linearities present after the caustic formation. Also, since 
it describes the evolution of a periodic, sinusoidal density pertur- 



The original pancake problem has an analytical solution {An-| 
|ninos & Norman 1994), which describes the collapse of a pres- 
sureless fluid up to caustic formation, happening at redshift z c (the 
moment of the first shell crossing). For baryonic collapse, it is valid 
as long as gas pressure is still negligible and can be used to set up 
the initial conditions. The Lagrangian positions and velocities are 
given by 
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Figure 10. Cosmological MHD pancake at z = with perpendicular magnetic field and dark matter for different initial field strengths B,,,,;. The test is run in a 
one-dimensional box with N = 4096 resolution for the hydrodynamical grid and 4 dark matter particles per grid cell. 
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To set up the test, we transform them to Eulerian coordinates: 

l+z c sin(£x,) 

x = x, - (46) 

1+z k 

This 'single pancake' test has been used by many authors to 
test cosmological hydrocodes, e.g. |Ryu et al.|fT993) , |Gheller et al.| 
|1996) , [Ricker et al.|p000) . We skip it here (it will appear again 



in section |3.1 1 and directly move on to the 'double pancake' test. 
It has been proposed by |Bryan et al.| C 1 995 ) and not been consid- 
ered by any other group ever since. In this test, a second wave with 
one fourth of the wavelength is superimposed on the original wave, 
utilizing the same formula {45}. The parameters of the two waves 
are 



A l = 64 Mpc/h 
Zcl = 1 



A 2 = 16 Mpc//? 
z c2 = 1.45 



and are evolved from z = 30 to z = 0. The initial baryon tempera- 
ture is set to T init = 13 K according to the formula given in |Anninos] 
|& Norman| ( [li96l >. 

The double pancake is a much more challenging test than the 
single pancake, not only because it introduces stronger shocks, but 
also because the superimposed additional wave leads to much finer 
features which are harder to resolve by the code, especially the 
temperature peaks. The ratio of thermal to kinetic energy density 
ps/^pv 1 in this test covers an extremely wide range between 1(T 9 
and 10 s , making it a very stringent test on the correct implementa- 
tion of the dual energy formalism. 



The numerical results are shown in Figure [8] While the low- 
resolution run with N = 64 fails to recover all features of the so- 
lution (the structure of the central density peak, the peak separa- 
tion in the temperature), they are present in higher resolutions. The 
high-resolution run with N = 1024 impressively recovers the so- 
lution of |Bryan et al.| ( |T995| >, and due to the dual energy method, 
the sharp side peaks in the temperature are resolved extremely well. 
Also, in the extremely cold regions outside of the peaks the temper- 
ature is tracked correctly without any oscillations or other artefacts. 
We could not reproduce this result without our dual energy imple- 
mentation or with other cosmological codes publically available. 
We acknowledge that the gravitational solver and the cosmological 
expansion (through supercomoving coordinates) are implemented 
correctly and that our variation of the dual energy formalism effec- 
tively improves tracking of the temperature without having a dy- 
namical effect on the density or velocity of the gas. 

We took a closer look at how well the solution of this test 
converges. For this, we ran the exact same test as described above, 
with different resolutions from A' = 8 to A' = 512. As long as the 
behaviour is linear, that is, well before caustic formation, we can 
define the relative A t error norm of a quantity q as: 

H N^f \q ref \ 

The left side of figure|9]shows this error for the density and velocity 
as a function of resolution at z = 20. We took the analytical solution 
(see above) as reference. Before the calculation of the error, the 
velocity was shifted by a constant, so that it does not approach zero. 
We find a constant convergence rate around N~ lA for the whole 
resolution range. 

For the final solution at z = 0, we are far in the non-linear 
regime, and the solution features strong discontinuities and even 
singularities in the density. If we are to make a similar study here, 
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we have to redefine what we take as the error. The analytical so- 
lution is not valid at this point, so we compare against a high- 
resolution run (N = 1024), which is binned down accordingly, and 
use the relative self-convergence error: 



I* _ 91024I 



N rnaxd^l, \q l0 2*\) 



(48) 



The denominator is chosen this way to force all terms to be between 
and 1. It leads to more meaningful results, because the differences 
can span over several orders of magnitude due to the strong discon- 
tinuities present. This error is plotted again for density and velocity 
at z = (the right side of figure|9]l. Because of the non-linearity of 
the system, the solution converges much slower at first and reaches 
N~ [ only at high resolutions. There is a minimum of resolution re- 
quired to get the shape of the solution right (around N = 128), 
for lower resolutions there are features missing. This is especially 
the case for the velocity distribution with its pronounced minima 
and maxima, producing the kink in the convergence rate between 
N = 64 and N = 128, and after that the convergence improves. 

When comparing this performance, it turns out that even in 
|Bryan et al. ( 1995), where the same test is run with a third-order 
accurate piece-wise parabolic (PPM) code (while our scheme is 
second-order piece-wise linear), the density distribution converges 
as A/~' 5 in the linear regime; and for z = 0, it does not get better 
than AT 1 either. The PPM code of Ricker et al^2000|l reaches only 



A' at z = 7 for the single pancake. In this context, we can safely 
state that our code performs adequately well. Being a second-order 
scheme, the MHD solver requires less computational steps and is 
faster than higher order methods, while attaining comparable accu- 
racy in the nonlinear regime of structure formation. 



3 COSMOLOGICAL MHD SIMULATIONS 

Having tested the functionality of the AMIGA code, we now move 
on and combine the A/-body solver, MHD and gravity to perform 
full cosmological MHD simulations. The aim of this section is to 
quantify the impact that the introduction of initial large-scale mag- 
netic fields has on simulations of the evolution of dark matter and 
baryons in a cosmological context. 



3.1 MHD pancake 

Before running simulations with realistic cosmological initial con- 
ditions, we use the Zel'dovich pancake collapse model in the sense 
of a single-mode analysis to get an idea of what effects to expect 
from the presence of a magnetic field. We use a wave in x-direction 
with z c = 1 and A = 64 Mpc/h, but with a few modifications. First, 
we also want to study the dark matter component. So we include 
dark matter particles and baryons simultaneously in the simula- 
tion, using a baryon fraction fi, = 0.165 (like later in the full sim- 
ulations). Both follow the same density and velocity distribution 
initially. Then, we apply a perpendicular, constant initial magnetic 
field Bj„i, pointing in (/-direction. This whole system is then evolved 
in a flat EdS universe (Q = 1, C1 A = 0) until z = 0. We repeat the 
same simulation for a wide range of different initial magnetic fields 
from B jnil = 10" 11 G, where the magnetic terms are neglectably 
small and dynamically unimportant, up to 6„„, = 2 • 10~ 6 G, where 
the magnetic field accounts for several percent of the total energy 
density of the gas. 

Figure [T0| shows the density and temperature of the baryons, 
the distribution of the dark matter particles and the magnetic field 




B Wt [G] atz = 30 

Figure 11. Relative deviation (Ai norm) of the final baryon and dark mat- 
ter density distribution from the B;,„, = case depending on B,,„, in the 
numerical MHD pancake solution. 



strength at z = for different runs. Initial fields up to up to 
0.05 /jG do not have any significant effect on the density pro- 
files (and the other quantities) and the field strength just follows 
the density profile of the baryons. Higher fields, however, induce 
large changes at the shock and post-shock regions, slowing down 
the baryon collapse and smearing out the baryon density profile. 
High field strengths effectively prevent the build-up of sharp, high- 
temperature baryon peaks. Although the situation is not directly 
transferrable to a full 3D simulation, this study gives us a hint on 
the general behaviour of density peak formation under the influence 
of a magnetic field. 

The dark matter distribution is generally much less affected 
than the baryon component - it does not interact directly with the 
magnetic field, but only indirectly through the gravitational force of 
the baryons. Since the dark matter particles are collisionless, they 
do not clump all together in the centre, but pass through each other 
(this happens exactly at z = z c ) and form side peaks. The shape of 
the central dark matter peak gets somewhat distorted if the baryons 
clump differently due to the field, but the height and position of the 
side peaks that form after z c lie outside the shocked regions and are 
practically not affected. 

For a more quantitative view on these effects, we took the non- 
magnetic numerical pancake solution at z - as a reference and 
calculated the deviation created by an initial magnetic field as the 
relative Ai "error" norm l |47| l of baryon and DM density distribu- 
tions over the interval of interest x e [0.43,0.57], a quantity very 
sensitive to numerical deviations (Figure [TT). Expectedly, the de- 
viation rises with higher initial fields, and the dark matter is less 
affected than the baryons. It should be noted that there is a cer- 
tain level of numerical noise: If one changes a code parameter, like 
the CFL number, initial timestep, or dual energy parameters (within 
reasonable values of course), the relative deviation will be typically 
around A t p as 10~ 4 . Smaller deviations can therefore be considered 
statistically insignificant. We can see that, in order to have a sig- 
nificant effect of, say, = 1%, the energy density of the initial 
magnetic field at z = 30 has to be at least around 10~ 7 G. 
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Figure 12. Comparison of cosmological MHD simulations at z = with low and high initial magnetic fields. Shown is the baryon density (top) and magnetic 
energy density (bottom) in code units. The maps are cuts through the simulation box in the x, y plane at half box depth. 

Unfortunately, the predictions of such models involving string the- 
ory and particle physics are presently highly parameter dependent 
and rather inconclusive. Certain models can lead to fields as large as 
1 nG, while others predict fields that are many orders of magnitude 
smaller ( Subramanian 2008 1. 

While at the present there is no working theory to estimate 
the possible strength of a primordial magnetic field, it is feasible 
to give some constraints. If a magnetic field coherent at cosmolog- 
ical scales was present in the early universe, it should have left its 
imprint on the linear polarization of the CMB by Faraday rotation. 
Based on that, Kahniashvili et al. (2009 ) derived an upper limit for 
a primordial magnetic field based on the CMB polarization power 
spectrum from the WMAP 5-year data. They find that at a scale of 
100 Mpc, the field amplitude must have been smaller than 0.7 nG. 
At smaller scales of 1 Mpc, the upper limit may be as high as 30 
nG, depending on the assumed power spectrum. 

Another way to constrain the primordial magnetic field stems 
from Big Bang Nucleosynthesis (BBN). The presence of a mag- 
netic field during BBN would have changed the nuclear reaction 
rates, thus resulting in an altered abundance of lighter elements 



3.2 Cosmological MHD simulations 

3.2.1 The initial magnetic field 

In order to conduct simulations of the cosmic structure forma- 
tion that take into account the primordial magnetic field, one must 
choose appropriate initial conditions. Yet, the possible strength^ 
shape or origin of a primordial field remains unclear at the mo- 
ment. Theories on this subject suggest that a cosmological large- 
scale magnetic field was already present before recombination (z ~ 
1000). Such a primordial magnetic field could have been produced 
during inflation (Turner & Widrow 1988; Gasperini 2006), much 
like the primordial density fluctuations that led to structure forma- 
tion, or during subsequent phase transitions (Gopal & Sethi 2005 1. 

4 When talking about primordial field strengths at earlier times (higher 
redshifts), we mean the comoving magnetic field strength B cmmv j„g = 
cT^Bproper, that is, the strength such a field would have when extrapolated 
to the present-day scale factor. It is convenient because it allows for the 
magnetic field strengths from different epochs to be compared directly. The 
relation B(z) = Bo/a 2 follows from magnetic flux conservation. 
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Simulation 


Initial comoving magnetic 


Initial physical ma; 


;netic 




field strength B,„;, at z = 30 


field strength at z = 


30 


BO 


OG 


OG 




Bl 


5.79 -10- 10 G 


5.56 -10- 7 G 




B2 


5.79 -10- 9 G 


5.56 lO^G 




B3 


5.79 -10- 8 G 


5.56 -10- 5 G 





Table 1. Identifiers and initial magnetic field strengths for the MHD simu- 
lations used in this paper 



like 3 He, 4 He, 7 Li. To be compatible with the current agreement be- 
tween BBN theory and element abundancy observations, a primor- 
dial field must be smaller than some critical value. First constraints 
derived in that manner were pretty high, up to 1 fiG jKernan et al.| 
1996); later, Grasso & Rubinstein (2001 1 deduced a more realistic 
value of 1 nG for Mpc-scale fields with the help of some additional 
assumptions. 

For the simulations presented in the following section, we as- 
sume a primordial field already present before the starting time of 
the simulation. We further assume it to be constant and homoge- 
neous in the whole simulation box, since the focus lies on how 
structure formation is affected by magnetic fields on scales larger 
than individual structures. It could be argued that a homogeneous 
primordial field pointing in one direction contradicts the assump- 
tion of an isotropic universe by creating a direction of preference; 
but the anisotropy created by such initial conditions has no impact 
on a statistical analysis of the baryon evolution, because the angle 
between the field vector and the baryon flow, a crucial quantity for 
the magnetic force on the baryons, is still randomly distributed, and 
the magnetic pressure does not depend on the field direction at all. 
It is also worth noticing that ideal MHD predicts the magnetic field 
lines to follow the baryon distribution. In fact, gravitationally col- 
lapsing baryonic structures would completely reshape the magnetic 
field distribution up to the point that any information on the original 
shape of the primordial field would be lost ( Dolag et al. 2 002] ), so 
the initial field shape should not significantly influence the results. 



3.2.2 Overview and initial conditions 

For this study, we carried out a set of cosmological 3D simulations 
with varying primordial field strengths S,„„. The simulations model 
a universe containing baryons and dark matter particles in a three- 
dimensional 64 Mpc//z box with periodic boundary conditions, us- 
ing a baryon fraction of f b = 0.165. The initial conditions used 
for all of the simulation runs were created from an initial CDM 
power spectrum corresponding to the WMAP 5 parameters (Ko- 
|matsu et al.|2009) : Q = 0.273 and Q. A = 0.726 with the PMCODE 
IC package (Klypin & Holtzman 1997). The initial density distri- 
bution of the baryons follows the dark matter. The initial field is 
set to a constant magnetic field in {/-direction: B,,„, = (0, 6,,„,,0). 
Apart from the different initial magnetic field strengths, the initial 
conditions are identical for all the runs. The simulations were run 
from the chosen starting redshift z,„„ = 30 until z = with the full 
MHD version of the AMIGA code on a regular 256 3 grid utilizing 
OpenMP statements for parallelization. 

The initial magnetic field values are summarized in Table 1. 
While the lowest initial field (Bl) is compatible with current theo- 
retical and observational estimates for a large-scale field, the high- 
est initial field (B3) is significantly higher than all current upper 
limits. The energy density of the B 1 field at z = 30 is equivalent 
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Figure 13. Power spectrum of baryons at z = for MHD runs with different 
initial magnetic fields relative to the power spectrum without magnetic field. 
For sufficient field strengths, the power spectrum of the baryon distribution 
shows a suppression of finer structures (the dark matter power spectrum 
stays unchanged). 





logp 



Figure 14. Contour plot of the volume fraction with given baryon temper- 
ature and density at z = without and with a strong initial magnetic field. 



to the kinetic energy density of a gas with the average gas density 
at that redshift moving at 0.4 km/s. This is clearly too small to be 
dynamically important at any stage of the simulation. The B3 field, 
on the other hand, is 10 2 times stronger and has 10 4 times more 
energy density, so we should expect an effect due to its presence. 



3.2.3 Magnetic field influence on the large-scale distribution of 
baryons 

Figure [T2| shows maps of the baryon density and magnetic energy 
density at the final redshift z = of our simulations, with the weak- 
est and strongest initial magnetic field that was simulated. Even for 
the very high magnetic field of the B3 run, the baryon distribution 
does not change significantly; a closer look reveals that the distri- 
bution becomes slightly smeared out, featuring less fine structure. 
The weaker initial fields do not have an influence at all. The mag- 
netic field shows the expected behaviour: in ideal MHD, it is frozen 
into the gas motion and largely follows the gas distribution. 

In the following, we want to analyse quantitatively whether 
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the presence of magnetic fields in these simulations alters the result 
for the other constituents. 

The smearing of structural features by magnetic fields can be 
quantified by looking at the power spectrum at z = 0. In figure [T3] 
we plot the power spectrum Pb ar (k) of the baryons relative to the 
power spectrum of the non-magnetic run. It can be seen that the 
additional magnetic pressure and the altering of baryon flows by 
magnetic tension leads to a characteristic suppression of structure 
formation at finer scales. This is, however, a very small effect that 
becomes important only at smaller scales and high magnetic fields. 
The high initial field in the B3 run lowers the power spectrum by 
10 % on a scale of 1 Mpc/h, the weaker field in the B2 model only 
has a (very slight) effect below 1 Mpc/h, and the more realistic 
field of the B 1 model has no effect whatsoever. On scales above a 
few Mpc//7, no influence can be seen even with the strongest field 
in our simulation runs. On scales well below 1 Mpc//i, additional 
processes inside the individual collapsing structures become im- 
portant, especially in their core regions. There, the magnetic field 
is influenced by cooling flows, turbulence and field tangling, which 
is not resolved by the large-scale runs conducted here. They can be 
better addressed by finer simulations of individual objects like the 
study of Dubois & Teyssier (20081 on a single magnetized galaxy 
cluster. 

Figure [14] shows the number of cells with a certain value of 
baryon density p and temperature T within the simulated box, again 
for the strongest initial field and without a field. While the lat- 
ter shows the characteristic shape known from other cosmological 
codes (see e.g. |Ryu et al.|19 93 for similar figures), a strong mag- 
netic field leads to an additional maximum at the bottom, where re- 
gions with very cool gas are located. We can compare this directly 
to the result for a single Zel'dovich wave in Figure[l0] where a sim- 
ilar effect occurs: a broadening and smearing of the density profile 
and the formation of a cool region behind the shock front. These 
results are in good agreement with Gazz ola~et al.| (|2007), where 
the same smoothing of the mass distribution with shallower den- 
sity profiles and "washed out" finer density clumps can be seen, al- 
though they add the magnetic field simply as an additional isotropic 
pressure term instead of a proper MHD treatment. In any case, for 
magnetic fields of order ~ 1 nG and below, no effect whatsoever 
can be seen on the scales resolved by our simulation. 

To summarize, on scales of ~ 1 Mpc and above, only primor- 
dial fields significantly higher than of order ~ nG have a noticeable 
effect on the baryon dynamics and gas distribution, which can be 
safely stated to be outside the upper theoretical and observational 
limits. 



4 SUMMARY AND CONCLUSIONS 

In this paper, we present the new numerical C code AMIGA de- 
signed to perform cosmological magnetohydrodynamic simula- 
tions. It contains the powerful and memory-efficient AMR A'-body 
code from its predecessor MLAPM Knebe et al. (2001 1, as well as a 
newly developed Eulerian grid-based MHD solver based on Ziegler 
(20041 and Ziegler (2005). The new code allows to simulate dark 
matter, baryon physics and magnetic fields in a self-consistent way 
inside a full cosmological framework. To facilitate the numeri- 
cal solution of cosmological MHD equations, the code is working 
with supercomoving coordinates, a transformation that greatly sim- 
plifies the equations, while preserving the fully cosmological set- 
ting. There are implemented techniques to properly resolve strong 
Shockwaves and supersonic flows in the baryon component, and to 



ensure the important condition of a divergence-free magnetic field 
down to machine precision. By conducting a series of test problems 
we acknowledge the high accuracy of this new code. 

As a first application of the new code, we present simulations 
of the cosmic structure formation with primordial magnetic fields. 
Such large-scale magnetic fields, possibly of cosmological origin, 
can be expected from different theoretical models and observational 
evidence of magnetic fields inside galaxy clusters. We want to ad- 
dress the question whether they could be a relevant factor for the 
large-scale dynamics in cosmological simulations. 

The simulations carried out with AMIGA model a ACDM uni- 
verse with the WMAP-5 cosmology in a comoving 64 Mpc//i com- 
putational volume, and its evolution from redshift Zinit = 30 to 
z = 0. The applied primordial field strengths range from about 
0.5 nG (a likely value from current constraints) to about 50 nG, 
which is significantly higher than current theoretical and observa- 
tional constraints for magnetic fields on such large scales. The anal- 
ysis of the simulations reveals that only in this last case, a large- 
scale magnetic field has a statistically significant influence on the 
baryon dynamics. Then, the magnetic pressure and tension leads to 
a suppression of baryonic small-scale structure and smears out den- 
sity peaks, visible in the baryonic power spectrum. However, even 
the highest simulated initial field has no noticeable effect on scales 
above a few Mpc//i. We can therefore conclude that, since current 
theoretical and observational constraints predict a large-scale field 
not much stronger than ~ 1 nG, at least outside of the core regions 
of gravitationally collapsed structures it cannot have a significance 
for the baryonic component during large-scale structure formation, 
neither on the power spectrum nor on the actual distribution. 

Even though our simulations do not have the required resolu- 
tion to study the (internal) properties of individual objects, we nev- 
ertheless like to close with a brief discussion of our findings in that 
direction. We observed (though not explicitly presented here) that 
the mass function of collapsed structures remains unaffected even 
for magnetic fields as large as the ones in model B3. Furthermore, 
the shape of (dark matter) haloes also appeared unaltered when in- 
creasing the strength of the primordial magnetic field. And for the 
baryon fraction - for which |Gazzola et al. (2007) have shown a 
dependence on the magnetic field strength - our own results are 
unfortunately affected by resolution effects: while stronger mag- 
netic fields lead to a depletion of baryons in smaller mass objects 
(cf. Figure 8 in Gazzola et al. ( 2007jl), the same is caused by a lack 
of resolution in cosmological codes jCrain et al.| 2007 ; Rud d et aT] 
2008); hence, we observe this effect but attribute it to our resolu- 
tion. Further studies and more refined simulations in this direction 
are necessary to clarify this subject in greater detail. 
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APPENDIX: DERIVATION OF THE SUPERCOMOVING MHD EQUATIONS 



As the starting point for the supercomoving transformation we take the full set of equations describing dark matter, baryons and magnetic 
fields in the ordinary non-comoving frame: 

dr dm 
dt 
dv dm 
dt 

A(p = ArcGpm, (Ale) 
dp 
dt 

d(pv) _ / B :: \ . 1 



+ V • (pv) = 



(Ala) 



(Alb) 



(Aid) 



+ V ■ 



+ V ■ 



dt 
d(pE) 

dt 
dB 

— + V x (-v x B) 

dt 

V B = 

dS 



v pE + p + 



B 1 



BB 



- P v<p 



2pJ p 



B(v B) 



-pv ■ (V0) 



dt 



+ V-(5u) = 



(Ale) 

(Alf) 

(Alg) 
(Alh) 

(Ali) 



Here, equations {ATa]l and jATb} are the equations of motion for the dark matter (dm) particles; equation (JaTcJ is Poisson's equation, where 
p w , is the total (dm + baryons) density; and the equations l |Aldfr - l |Alh[ > are the equations of ideal MHD. For the dual energy formalism, we 
also need to transform the equation {Ali} describing the modified entropy. Below, we will apply the supercomoving transformation to each of 
these equations individually and construct a new set of supercomoving equations, using the definitions |TJ and l|2j. The new supercomoving 
quantities and derivatives will be denoted by a subscript x. Throughout this appendix, an overdot denotes the temporal derivative with respect 
to the proper, non-supercomoving time t. 



4.1 Dark matter particle equations of motion 

We define the supercomoving velocity v x j„, as the derivative of xj,,, with respect to the supercomoving time df v = dt/a 2 



dx 



dill 



' vx dm 



From this follows the relation between physical and supercomoving velocity: 



^dm ^x,dm @%dm 

a 



(A2) 



(A3) 



Qx,dm Q^dm Q^dm '■' 

It follows also from the definition that the spatial derivatives change to 

V_ v = aV ; A_ v = a 2 A . 

The goal is to obtain an equation of motion analogous to equation [Alb] for the supercomoving dark matter velocities v x ^m ■ Let us consider 
the supercomoving acceleration: 

dVdm_ _ d l dx dm \ 
dt x dt x \ dt x J 

Now we replace the supercomoving time and position with the physical time and position, and rewrite the result to obtain a relation with the 
supercomoving gravitational force: 



dv x . 



d/, 



- a —(V x dm) 

dt 



: a 2 — (av dm - ar dm ) 

- 1 dVd„, 
n 

dt 

dvd„, 

i 

dt 

dVdm 
I 

dt 



dt 

= a' | a — — + av im - ar dm - av dr . 



-V,cA + d"[-^ x \^X 2 dm 
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= -v„ 

= -Vx<l>x 



1 .. 2 

> + -aax rh 



Here, we used equation ( |Alb[ > and the definition of the comoving potential (f> x . Now we see that the supercomoving equations of motion are 
formally identical to their proper physical counterparts, although the quantities are defined differently: 

ted,,, , . . . 

—. — = i>*,dm (A4a) 

-jf* = -V,0.v (A4b) 

This is the main advantage of supercomoving coordinates over the comoving coordinates, which explicitly include additional factors depend- 
ing on a. We will see that the other equations behave in a similar way. 

4.2 Poisson's equation 

Poisson's equation determines the potential if> of the system, and as such it is the only equation where cosmology enters explicitly. Although 
equation < |Alc[ l describes the gravitational potential in an ordinary physical setting, in a cosmological framework we also have to consider 
the cosmological constant A. This is realized by adding a A term that has the dimension of a density, and then using this "effective density" 
in Poisson's equation. 

Let us consider the second Friedmann equation, which relates the average total mass density p m and the cosmological constant to the 
accerelation of the cosmic expansion: 

a AnG I 3»\ Ac 2 

a=-—[^ + i) + - (A5) 
Dark matter is by definition pressureless (p = 0), and the pressure of the small baryonic component can be neglected. Then we can write 

- = - — r~ (fitot ~ Pa) (Ad) 
a 3 

with pa = -he 1 1 AnG. Then, if the cosmological constant is not zero, Poisson's equation effectively becomes 

Acf> = 47rG(p tol -p A ) . (A7) 
We formulate the left-hand side in terms of the supercomoving gravitational potential: 

1 (<f> x 1 2 

A(4 = —A, - -aax 
a 1 \a l 2 

1 d 2 

= ~zh x <j> x - —A x x 
a La 

1 i , 



CT a 



Using again the second Friedmann equation|A6| we get 



= -jA,0, - 3 
a 4 



-AnG 

— z — (Pwt- Pa) 



Equating this with the right-hand side of equation |A7 



— ■ A r f r - AnG(p m - p A ) = -AnG(p m - p A ) . 
a* 

The A term cancels, and we arrive at the supercomoving Poisson's equation: 

A. r 0.v = AnGa A (p m - p m ) 

= AnGa(p xm - p xm ) . (A8) 

The supercomoving version of Poisson's equation looks slightly different than the non-cosmological one: the density contrast enters instead 
of the total density, because the supercomoving potential is responsible for peculiar motions due to density fluctuations, while the total density 
governs the overall expansion. 

4.3 Baryon mass density 

Now we will transform the conservation law for the baryon density p, 

dp 
at 
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First, we replace proper time, density and velocity with their comoving counterparts. We begin by replacing the partial time derivative d/dt 
at constant position r with the one at constant comoving position x: 



P 



p - -x ■ V v p 

a 



d_ 
dt 
1 

a 2 dt x \a 3 
-3a 

a 



1 d I 1 \ a „ 

3 A ~ 7a x ' VxPx 



1 dp x a 



a 3 dt x ar 



Evaluating the flux term, using v = v x /a + ax: 



1/1 a 
v(po) = -V,| -^p x V x + —p x x 



^x(p,v x ) + -rV x (p x x) 
^ a 4 



a \ar 
1 

1 t, , x a ,„ ^ 3a 
= -rYrfaA) + ^x(V v p r ) + —p x 
a? a 4 a 4 

Putting these two expressions together, four of the six terms cancel, leaving only 

dp. 



+ VAp*v x ) = 



(A9) 



4.4 Baryon momentum density 

The conservation law for the baryon momentum can be written as 

Bipv) 



dt 



+ V • (pvv + A) = -p V(p 



(A10) 



where we used the abbreviation A = (p + 0) / - j)BB. From the definitions of the supercomoving pressure p x and magnetic field B x one 
immediately sees that 



A, = [p x + ^j/- jB x B x = a 5 A 



We decompose the first term of the momentum equation into two parts 

dipv) 



dt 



dv 


dp 







The first part equals 



dv 

'ft 



1 / dv 

a~M P Jt 



-x(V x u) 



Using the abbreviation K = c jx(V x ii) , it evaluates to 

dv 

P di 



1 / dv 

a 3P \ P dt 



- K 



1 



x \ a 



-v x + ax - K 



while we write the second part as (equation |Alb[ l: 



dt 



1 dv x 

a 3 dt x 



-uV • (pv) 



The second term on the left-hand side of equation l[AT0} transforms as follows: 

V(pi>u + A) - uV • (pv) + pv • V« + VA 



= vW ■ (pv) + —p x 

a s 



vV ■ (pv) + —p x 

a s 



vW ■ (pv) + —p x 

a 3 



1 \ 1 / 1 \ 

-v x + ax • - V v —V x + ax + K 
a J a \a j 

\(v x -V x )v x + ^-(v x -V x )x + K 

a 1 a 1 



a" 



-^(P x ■V x )v x + ^-v x + K 
a 1 or 



1 



a b 
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Combining all terms from the left-hand side of equation (JaToJ gives 



d(pv) 



Bt 



+ V(pvv + A) = p 



dv 
7k 



dp 

'a? 



+ V(pvv + A) 



1 



1 do x 1 



a 3 dt x 



+ — (v x -V x )v x + iix 



a 6 



When comparing this to the right-hand side of equation dATOl, 



ar \a- 



-1 



1 



-Px(.^x<t>x) + ^rP.vV., -ax 



we notice that the {ax) term cancels, leaving 



Px^T+Px(l>x ■ Vx)v* + Vx A x = ~Px^x4>x 

at x 



d(p x v x ) 

dt x 
d(p x v x ) 

dt x 



+ V,(/),B,B, + A x ) = -p x V x <p x 
+ V. t 



j bV\ i 

P X V X V X + \Px+7r\ I ~ ~ B * B * 
2pJ n 



-Px^xl 



(All) 



4.5 Induction equation 

Before transforming the total energy density equation in the next subsection, the supercomoving induction equation has to be derived as it 
will be needed for it. We start from the equation 



dB 

7k 



subject to the condition 



+ V x (-» x B) = 



V B = 



+ V x (-0 X B) = 



ix 1 
- -x(V x -B) + -V v x 
a a 



- 1 -v x + ax \ x B 







The divergence-free condition turns out to be useful as it causes several terms to vanish. First we substitute the temporal and spatial deriva- 
tives: 

dB 

dt 
dB 

^ ~dl 

<=> 4f - 4v,x(bxB)--V,x(xxB) = 

a 1 at x a 1 a 

Again, the divergence-free condition allows us to simplify: 

V v x (x x B) = (B V)x - B(V ■ x) = B - 3B = -IB 

and therefore 

1 <>B 1 V t x(uxB) + 2-B = 
a 



a 2 dt x a 2 
dB 



dt, 



V x x(vxB) + laaB = 



Now we substitute the supercomoving magnetic field B x = a 5/2 B: 



a -5/2?Bx + Bx ±i-S,2\ _ a -5,2 Vjc x ( x g x) + lha -V2 Bx = Q 

ot x dt x v ' 

5/2 ^£ -5/2 ' 



a~ 3// V v x(v x XB x )- -aa- 3/2 B x = 
at, 2 



We define the supercomoving Hubble constant 



1 da 

li := - — = aa 
a at x 
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With this notation, we have: 

+ V, x (-», x B x ) = \<HB X (A 1 2) 

We defined the frame of reference such that it is comoving with the magnetic energy density, and not with the magnetic field strength. 
This is the reason why a magnetic Hubble drag term must appear at the right-hand side of the supercomoving induction equation and it is not 
formally identical to the non-comoving induction equation. However, the Hubble term only ensures that the magnetic field scales properly 
with a; it does not have any physical meaning. 



4.6 Total energy density 

Instead of directly transforming the total energy equation l |Alf) , we derive the supercomoving energy conservation law by putting together all 
the quantities we have so far. The easiest way is to first derive the energy conservation for the hydrodynamic case and then add the magnetic 
energy density and flux to the result. 

In the hydrodynamic case, pE = jpv 2 + ps. One immediately notices from the definitions of the supercomoving variables that the 
supercomoving total energy is 

1 2 

PxE x = -jPxV x + p x e x 

We start by calculating the temporal change of the kinetic energy ^p x v\ with the help of the already derived equations: 

d (1 



dv x 1 2 dp x 



ft, 1 2 



1 



1 



-v x V x (p x v x ) 



- p v (v x ■ V )v. - V (.'■, V /> 

i , 

; PxVx [-(»* • VjDj - ^V X V X ■ (p x V x ) - V x ■ ~ x p x -p x V x ■ (\ x <p x ) 

: -V.v • (ox-p^j - v x ■ V,/), -p,v, ■ (V x (p x ) 



With v x ■ p x = V v ■ (p x v x ) — p x V x ■ v x , we can write: 



d (1 



dt x 12 



/','•. +V, 



1 ~pA + p. <]•■>■ 



■ Px^x ■ V x - PxVx ■ (Vx<Px) 



(A 13) 



Next, we need an equation for the thermal energy e x . In proper coordinates, such an equation exists (e.g. |Bryan et al.|1995| . In the case of a 
monoatomic ideal gas (7 = 5/3), which will be assumed from here on, it reads: 



ds 1 

— + V ■ \E = p\ ■ V 

Ot p 



(A14) 

By plugging in the definitions of the supercomoving variables it easily proves that the same equation holds in supercomoving coordinates: 

(A 15) 



de x 1 „ 

— + v x ■ V x e x = p x S x ■ v x 

dtx Px 



We can rewrite that as 



„ de x _ 

P"x ■ v x = ~Px^ PVx ■ V_ v e v 

dtx 



S(pxEx) 

dtx 



- V x [v x (p x e x )] 



and plug it into equation l |A13[ l, yielding 



dt x 



(p x E x ) + V, • [ifxEx + Px) v x ] = -PxVx ■ (V x <p x ) 



(A 16) 



This is the supercomoving total energy equation for the hydrodynamic case, again formally equivalent to the corresponding non-comoving 
equation. 

Now we can consider the magnetic energy B\j2p. Its temporal derivative is easily obtained from the supercomoving induction equation: 

dtx\2p-j p 1 dtx 



1 



V v x (v x x B x ) + ^HB X 



1 . 



-V^ x (v x x B x ) 



2p 
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V, 



(v x XB x )xB x 



2p 



I S 2 
-fl,(o. t -B,)-u, -i 



Adding this equation to < |A 1 6^ , we get the supercomoving total energy equation for the full MHD case: 



8t x 



(p x E x ) + V, 



(fi x E x + p x )v x - -B x (v x ■ B x ) 
ft 



in 



where now 



1 2 B \ 

PxE x = -p x v x +p x e x + — . 

4.7 Modified entropy 

This additional equation is needed to use the "S system" in the dual energy formalism. Transforming the first term: 



(A17) 



and the second term: 



Putting both together yields 



dS 


dS 








P 

X 




i 


8 




a 2 dt x 




= a 3 >~ 


10 



p--x-V,S 

a 

(a 3 ^S x )-^x-V x (a^- s S x ) 
dS x 



+ 'H(3y- 8)S x -<Hx-V x S x 



V ■ (Sv) = -V, 
a 



a 3y S S X \ —V x + ax 



dS x 
dt x 
8S X 
dt x 



a 3r - l0 V x -[S x v x +'HxS x ] 

+ "H(3y - 8)5. v + V v • (S x v x ) + 39iS x = 
+ V,.(5>*) = --W(3y-5) . 



(A 18) 



